% =====================================================================
% ====  "Risk-Sharing in Village Economies Revisited", by Bold and Broer

% ====   This programme choses the preference parameters (CRRA rho,
% discount factor beta) according to various criteria

% ====    Written by Tessa Bold and Tobias Broer

% ====    This version June 2021
% =====================================================================%
clear
clf


      
Momvec1data=NaN(18,1);
Momvec2data=NaN(18,1);
Tessa1Tobi2=1;
outlier=0;
analytic=0;
est_rho1=0;
unconditional=0;
village=1
COALDEV=[0 1];
% =====================================================================%
% ====    Data moments
% =====================================================================%

clear icrisatdata momentsdatabootstrap DCdata DYdata DCdataneg DYdatapos DCdataneg DYdataneg Cdata Ydata VCV1data VCV2data
clear RelCYVardata RelDCDYVardata betaDCDYdatadyneg betaDCDYdatadypos betaDCDYdata Momvec1datamat Momvec2datamat DCskewdata RelDCDYVardata Cdata Ydata DCdata DYdata
clear betaDCDYdata CVARdata YVARdata DCVARdata DYVARdata Cskewdata DCskewdata RelCYVardata RelDCDYVardata


    
if Tessa1Tobi2==1
    outputpath='C:\Users\tbold\Dropbox\Tessa and Tobi\JEEA Replication\Programmes\APPROXIMATE_EXACT';
    datapath='C:\Users\tbold\Dropbox\Tessa and Tobi\JEEA Replication\Data/Output';
    programpath='C:\Users\tbold\Dropbox\Tessa and Tobi\JEEA Replication\Programmes';
else
  
    outputpath='C:\Users\tbroer\Dropbox\Research\Current_Projects\tessa\JEEA Replication\Programmes\APPROXIMATE_EXACT';
    datapath='C:\Users\tbroer\Dropbox\Research\Current_Projects\tessa\JEEA Replication\Data/Output';
    programpath='C:\Users\tbroer\Dropbox\Research\Current_Projects\tessa\JEEA Replication\Programmes';  
end


datafile=['icrisatmoments' int2str(outlier)];

    
    cd(datapath)
    
    eval(['load ' datafile '.out'])
    eval(['datamoments=' datafile ';']);
    eval(['clear ' datafile])
    
    
    Nom=4;
    
    bootstrapvar=datamoments(1:Nom,4+(village-1)*Nom+1:4+(village)*Nom);
    analyticvar=datamoments(1:Nom,16+(village-1)*Nom+1:16+(village)*Nom);
    analyticmoments=datamoments(1:Nom,1:3);
    
    analyticrawmoments=datamoments(1:Nom,28+1:28+3);
    bootstraprawvar=datamoments(1:Nom,28+4+(village-1)*Nom+1:28+4+(village)*Nom);
    analyticrawvar=datamoments(1:Nom,28+16+(village-1)*Nom+1:28+16+(village)*Nom);
    %%%here you replace the exact standard errors (the variances are
    %%%bootstrapped in STATA, the coefficients are delta-method
    
    VCV2conddata=bootstrapvar;
    VCV2rawdata=bootstraprawvar;
    
    if analytic==1
        VCV2conddata(2,2)=analyticvar(2,2);
        VCV2conddata(4,4)=analyticvar(4,4);
        VCV2rawdata(2,2)=analyticrawvar(2,2);
        VCV2rawdata(4,4)=analyticrawvar(4,4);
    end
    
    
    
   
    Momvec2conddata=analyticmoments(:,village)';
    Momvec2rawdata=analyticrawmoments(:,village)';
    if unconditional==1
        VCV2data=VCV2rawdata;
        Momvec2data=Momvec2rawdata;
    else
             VCV2data=VCV2conddata;
        Momvec2data=Momvec2conddata;  
    end
    
    COUNT=0;
for coaldev=COALDEV
        for exact=1:2
        cd(outputpath)
            
        eval(['load COMP_Momentsest_repl' int2str(coaldev), int2str(exact) '.mat'])
        
        cd(programpath)
        vss_set=[2 3];
        
        betavec=[squeeze(Momentsest(2,exact,:,1,2))];
  
      
       for inddelta=1:length(betavec)

             for vss=vss_set
             Vardylog=Momentsest(22,exact,inddelta,1,vss);              
             Vardylog_cond=Momentsest(30,exact,inddelta,1,vss);
                            

                            Momentstable(1:3,inddelta,vss,coaldev+1,exact)=Momentsest(1:3,exact,inddelta,1,vss);
                            Momentstable(4:7,inddelta,vss,coaldev+1,exact)=[Momentsest([8 15],exact,inddelta,1,vss) ; Momentsest(9,exact,inddelta,1,vss)/Vardylog_cond(village) ;(Momentsest([16],exact,inddelta,1,vss)-Momentsest([17],exact,inddelta,1,vss))];
                            Momentstable(8,inddelta,vss,coaldev+1,exact)=(Momentsest(21,exact,inddelta,1,vss)-Momentsest(21,exact,inddelta,1,vss))./Momentsest(21,exact,inddelta,1,vss);%2*maxcerr*(j-1)/(size(moments,2)-1)/moments(28,j);
                            Momentstable(9,inddelta,vss,coaldev+1,exact)=2*Momentsest(5,exact,inddelta,1,vss)./Vardylog;
                                                
                            Momentstable(1,inddelta,vss,coaldev+1,exact)=Momentstable(1,inddelta,vss,coaldev+1,exact)+1;
                            momentsdev=Momentstable(4:7,inddelta,vss,coaldev+1,exact)- Momvec2data(1,1:4)';
                            %%%2,4, d:diagonal, 
                            
                            W4(inddelta,vss,coaldev+1,exact)=momentsdev'*inv(VCV2data(1:4,1:4))*momentsdev;
                            W4d(inddelta,vss,coaldev+1,exact)=momentsdev'*inv(diag(diag(VCV2data(1:4,1:4))))*momentsdev;
                            W4_eye(inddelta,vss,coaldev+1,exact)=momentsdev'*inv(eye(4))*momentsdev;

                                                    
                            
             end
       end
    %%Estimation                    
      WW=W4d(:,vss_set,coaldev+1,exact);
      criterion=reshape(WW,prod(size(WW)),1);
      [CC,II]=min(criterion);
      [II1,II2]=ind2sub(size(WW),II);
      solution_4(:,coaldev+1,exact)=[Momentstable(1:9,II1,vss_set(II2),coaldev+1,exact); CC];
       
      WW=W4_eye(:,vss_set,coaldev+1,exact);
      criterion=reshape(WW,prod(size(WW)),1);
      [CC,II]=min(criterion);
      [II1,II2]=ind2sub(size(WW),II);
      solution_4_eye(:,coaldev+1,exact)=[Momentstable(1:9,II1,vss_set(II2),coaldev+1,exact); CC];
   
   
    end                  
                      
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%inspection CD
Momentstable(4:7,:,2,2,1)
Momentstable(4:7,:,2,2,2)
Momentstable(4:7,:,3,2,1)
Momentstable(4:7,:,3,2,2)
%%%%inspection ID
Momentstable(4:7,:,2,1,1)
Momentstable(4:7,:,2,1,2)
Momentstable(4:7,:,3,1,1)
Momentstable(4:7,:,3,1,2)




%%%Inspection by hand, picture for paper down below.
%%%average variance, group of size 3, CD: exact and approximate, ID: exact
%%%and approximate solution
subplot(2,2,1)
plot(betavec,squeeze(Momentstable(4,:,2,2,2)),'-b')
hold on
plot(betavec,squeeze(Momentstable(4,:,2,2,1)),'--b')
plot(betavec,squeeze(Momentstable(4,:,2,1,2)),'-r')
plot(betavec,squeeze(Momentstable(4,:,2,1,1)),'--r')
%%%average variance group of size 4, CD: exact and approximate, ID: exact
%%%and approximate solution
%%%solution
subplot(2,2,2)
plot(betavec(3:end),squeeze(Momentstable(4,3:end,3,2,2)),'-b')
hold on
plot(betavec(3:end),squeeze(Momentstable(4,3:end,3,2,1)),'--b')
plot(betavec(3:end),squeeze(Momentstable(4,3:end,3,1,2)),'-r')
plot(betavec(3:end),squeeze(Momentstable(4,3:end,3,1,1)),'--r')

%%%average degree of risk-sharing, group of size 3, CD: exact and approximate, ID: exact
%%%and approximate solution
subplot(2,2,3)
plot(betavec,squeeze(Momentstable(5,:,2,2,2)),'-b')
hold on
plot(betavec,squeeze(Momentstable(5,:,2,2,1)),'--b')
plot(betavec,squeeze(Momentstable(5,:,2,1,2)),'-r')
plot(betavec,squeeze(Momentstable(5,:,2,1,1)),'--r')
%%%average degree of risk-sharing group of size 4, CD: exact and approximate, ID: exact
%%%and approximate solution
%%%solution
subplot(2,2,4)
plot(betavec(3:end),squeeze(Momentstable(5,3:end,3,2,2)),'-b')
hold on
plot(betavec(3:end),squeeze(Momentstable(5,3:end,3,2,1)),'--b')
plot(betavec(3:end),squeeze(Momentstable(5,3:end,3,1,2)),'-r')
plot(betavec(3:end),squeeze(Momentstable(5,3:end,3,1,1)),'--r')


%%%asymmetry moments, group size 3, exact, CD, ID
figure(2)
%%%asymmetry variance, group size 3, CD: exact and approximate, ID: exact
%%%and approximate
subplot(2,2,1)
plot(betavec,squeeze(Momentstable(6,:,2,2,2)),'-b')
hold on
plot(betavec,squeeze(Momentstable(6,:,2,2,1)),'--b')
plot(betavec,squeeze(Momentstable(6,:,2,1,2)),'-r')
plot(betavec,squeeze(Momentstable(6,:,2,1,1)),'--r')
%%%asymmetry variance, group of size 4, CD: exact and approximate, ID: exact
%%%and approximate solution
%%%solution
subplot(2,2,2)
plot(betavec(3:end),squeeze(Momentstable(6,3:end,3,2,2)),'-b')
hold on
plot(betavec(3:end),squeeze(Momentstable(6,3:end,3,2,1)),'--b')
plot(betavec(3:end),squeeze(Momentstable(6,3:end,3,1,2)),'-r')
plot(betavec(3:end),squeeze(Momentstable(6,3:end,3,1,1)),'--r')

%%%asymmetry degree of risk-sharing, group of size 3, CD: exact and approximate, ID: exact
%%%and approximate solution
subplot(2,2,3)
plot(betavec,squeeze(Momentstable(7,:,2,2,2)),'-b')
hold on
plot(betavec,squeeze(Momentstable(7,:,2,2,1)),'--b')
plot(betavec,squeeze(Momentstable(7,:,2,1,2)),'-r')
plot(betavec,squeeze(Momentstable(7,:,2,1,1)),'--r')
%%%asymmetry degree of risk-sharing group of size 4, CD: exact and approximate, ID: exact
%%%and approximate solution
%%%solution
subplot(2,2,4)
plot(betavec(3:end),squeeze(Momentstable(7,3:end,3,2,2)),'-b')
hold on
plot(betavec(3:end),squeeze(Momentstable(7,3:end,3,2,1)),'--b')
plot(betavec(3:end),squeeze(Momentstable(7,3:end,3,1,2)),'-r')
plot(betavec(3:end),squeeze(Momentstable(7,3:end,3,1,1)),'--r')


%%%%%this produces Figure 7 in Section A.5.1. 

X1=betavec';
YMatrix1=([Momentstable(4:7,1:end,2,2,2)-Momentstable(4:7,1:end,2,2,1)])'; %CD, GS=3, exact, approximate
YMatrix3=([Momentstable(4:7,3:end,3,2,2)-Momentstable(4:7,3:end,3,2,1)])'; %CD, GS=4, exact, approximate
YMatrix2=([Momentstable(4:7,:,2,1,2)-Momentstable(4:7,:,2,1,1)]); %ID, GS=3, exact, approximate
YMatrix4=([Momentstable(4:7,:,3,1,2)-Momentstable(4:7,:,3,1,1)]); %ID, GS=4, exact, approximate

exact_approximate_picture(X1, YMatrix1, YMatrix2, YMatrix3, YMatrix4)
print('exact_approximate_pic','-dpdf','-bestfit')



